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Mass distribution of highly flattened galaxies and modified Newtonian dynamics 


W. F. Kao 

Institute of Physics, Chiao Tung University, HsinChu, Taiwan 

Dynamics of spiral galaxies derived from a given surface mass density has been derived earlier in 
a classic paper. We try to transform the singular elliptic function in the integral into a compact 
integral with regular elliptic function. Solvable models are also considered as expansion basis for RC 
data. The result makes corresponding numerical evaluations easier and analytic analysis possible. It 
is applied to the study of the dynamics of Newtonian system and MOND as well. Careful treatment 
is shown to be important in dealing with the cut-off of the input data. 
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I. INTRODUCTION 

The rotation curve (RC) observations indicates that less than 10 % of the gravitational mass can be measured 
from the luminous part of spiral galaxies. This is the first evidence calling for the existence of un-known dark matter 
and dark energy. In the meantime, an alternative approach, Modified Newtonian dynamics (MOND) proposed by 
Milgrom Hi: has been shown to agree with many rotation curve observations [l4 ll5L Il7t 1771 I 20 L l2~lT| . 

Milgrom argues that dark matter is redundant in the approach of MOND. The missing part was, instead, proposed 
to be derived from the conjecture that gravitational field deviates from the Newtonian 1/r 2 form when the field 
strength g is weaker than a critical value go ~ 0.9 x 10~ 8 cm s -2 [2ll| . 

In MOND the gravitational field is related to the Newtonian gravitational field gw by the following relation: 


9 ' M —) = 9n 
9o 

with a function go considered as a modified inertial. Milgrom shows that the model with 

g 0 {x) = 


\/l + x 2 


(1) 


( 2 ) 


agrees with RC data of many spiral galaxies iimmoiiiimEii The alternative theory could be compatible with 
the spatial inhomogeneity of general relativity theory. Various approaches to derive the collective effect of MOND 
has been an active research interest recently, jltj, 1221 

Recently, it was shown that a simpler inertial function of the following form sm 


9 o(x) = J 


+ X 


( 3 ) 


fits better with RC data of Milky Way and NGC3198. Indeed, one can show that the MOND field strength g is, from 
Eq. (1), 


9 = 


V 9 n 2 + 4<7q<7jv + gN 


( 4 ) 


The model shown by Eq. (2) will be denoted as Milgrom model, while the model (3) will be denoted as FB model. 
We will study and compare the results of these two models in the paper. 

Accumulated evidences show that the theory of MOND is telling us a very important message. Either the Newtonian 
force laws do require modification in the weak field limit or the theory of MOND may just represent some collective 
effect of the cosmic dark matter awaiting for discovery. In both cases, the theory of MOND deserves more attention 
in order to reveal the complete underlying physics hinted by these successful fitting results. Known problem with 
momentum conservation has been resolved with alternative covariant theories[217t ilOL Ull ll3LII 81 1221 . 

It was also pointed out that the proposed inertial function g 0 could be functional of the whole N field and could 
be more complicate than the ones shown earlier mm Successful fitting with the rotation curves of many existing 
spiral galaxies indicates, however, that the most important physics probably has been revealed by these simple inertial 
functions shown in this paper. One probably should take these models more seriously in order to generate more clues 
to the final theory. 

If the theory of the MOND is the final theory of gravitation without any dark matter in a large scale system 
with some inertial function go(,x), one should be able to derive the precise mass distribution form the measured RC 
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provided that the distance D is known. Earlier on, the mass distribution derived from the rotation curve measurement 
can only be used to predict how much dark matter is required in order to secure the Newtonian force law. In the case 
of MOND, one should be able to plot a dynamical profile of the E function T(r) = M(r)/L(r ) with the measured 
luminosity function L(r) following the theory of MOND. The P function should then provide us useful information 
about the detailed distribution of stars with different luminous spectrum obeying the well-known L ~ M a relation. 

The dynamical profile T(r) can be treated as important information regarding the detailed distribution of stars with 
different mass-luminosity relations within each spiral galaxy[s|. Even one still does not know a reliable way to derive 
the inertial function no(x), there are successful models shown earlier as useful candidates for the theory of MOND. 
It would be interesting to investigate the dependence of these models with the M/L profiles. Hopefully, information 
from these comparisons will provide us clues to the final theory. 

Therefore, we will try first to derive the mass density for Milgrom and Famaey&Binney (FB) models in details. 
There are certain boundary constraints needed to be relaxed for the asymptotic flat rotation curve boundary condition 
which is treated differently in the original derivation of these formulae (23 • It is also important to compare the effects 
of exterior contribution between the Milgrom and FB models for a more precise test of the fitting application. 

We will first review briefly how to obtain the surface mass density /i(r) from a given Newtonian gravitational field 
gisr with the help of the elliptic function K (r) in section II. The integral involving the Bessel functions is derived in 
detailed for heuristical reasons in this section too. 

In addition, a series of integrable model in the case of Newtonian model, Milgrom model as well as the FB model 
for MOND will be presented in section III. These solvable models will be shown to be good expansion basis for the 
RC curve data for spiral galaxies. In practice, this expansion method will help us better understand the analytic 
properties of the spiral galaxies. 

We will also try to convert the formula shown in section II into a simpler form making numerical integration more 
accessible in section IV. The apparently singular elliptic function K(r) is also converted to combinations of regular 
elliptic function E(y) by properly managed integration-by-part. 

In section V, one derives the interior mass contribution /r(r < R) from the possibly unreliable data v(r > R) both 
in the cases of MOND and in the Newtonian dynamics. The singularity embedded in the useful formula is taken care 
of with great caution. Similarly, one tries to derive the formulae related g]y from a given /r(r) in section VI. One also 
presents a simple model of exterior mass density /i(r > R) in this section. The corresponding result in the theory of 
MOND is also presented in this paper. Finally, we draw some concluding remarks in section VII. 


II. NEWTONIAN DYNAMICS OF A HIGHLY FLATTENED GALAXY 


Given the surface mass density /i(r) of a flattened spiral 
convert /z(r) to fi(k) in fc-space via the following equations 


galaxy, one can perform the Fourier-Bessel transform to 

0 


H(r) 

/i(k) 



kdk g,(k)Jo(kr), 


rdr /i(r) Jo(fcr) 


(5) 

( 6 ) 


with J m (x) the Bessel functions. Note that the closure relation 



kdkJ m (kx)J m {kx') = —5(x — x') 
x 


(7) 


can be used to convert Eq. © to Eq. © and vice versa with a similar formula integrating over dr. The Green 
function of the equation 


X7 2 G(x) = -4irG6{r)6(z) 


( 8 ) 


can be read off from the identity 


Vr 2 + 2 2 


dk exp[— k\z\] Jo(fcr). 


Therefore, the Newtonian potential (f>N can be shown to be [Hill 

pOO 

/ dk fi(k)Jo(kr) exp[—k\z\] 


(9) 


0N(r,z) = 27rG 


o 


( 10 ) 
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with a given surface mass density g{r). It follows that 

poo 

<7iv(r) = -~d r (f>N(r, z = 0) = 2ttG / kdk /i(k)Ji(kr). 

Jo 


( 11 ) 


The Newtonian gravitational field ^(r) = v n( r ) / r can be readily derived with a given function of rotation velocity 
Vi v(r). Here one has used the recurrence relation J\(x) = —J' 0 (x) in deriving Eq. (Hill . Furthermore, from the closure 
relation 0 of Bessel function, one can show that the function ffAr(r) satisfies the following conversion equation 


Therefore, one has 


pOO poo 

9n{t) = / kdk r'dr' g N {r')Ji(kr)Ji(kr'). 

J 0 Jo 

J rdr g N (r) Ji(fcr). 


Hence, with a given surface mass density /i(r) for a flattened spiral galaxy, one can show that 

1 


POO POO 

g{r) = —tG J kdk J dr ' r '9N(r')J 0 (kr)J 1 (kr'). 


( 12 ) 


(13) 


(14) 


Assuming that linv^o r ffiv( r ) —* 0 and lim^oo rg^(r) < oo hold as the boundary conditions, one can perform an 
integration-by-part and show that 


i p OO p OO 

^( r ) = 27tG J dk j dr ' d r'[r'g N (r')\J 0 (kr)J 0 {kr') (15) 

with the help of the asymptotic property of Bessel function: J m (r —> oo) —> 0. Here we have used the recurrence 
relation J\{x) = —J' 0 (x) in deriving above equation. 

Note that the boundary terms can be eliminated under the prescribed boundary conditions. In fact, the limit 

lim r _ rgjsr^r'j < cxd is slightly different from the original boundary conditions given in l23l . The difference is aimed 

to make the system consistent with the flatten RC measurement which implies that linv^oo r^Ar(r) = v^(r —» oo) —> 
constant < oo. In summary, one only needs to modify the asymptotic boundary condition in order to eliminate the 
surface term during the process of integration-by-part. We have relaxed this boundary condition to accommodate 
any system with a flat rotation curve. 

One can further define 



dkJo(kr)Jo(kr') 


(16) 


and write the function g(r) as 


mM = 


2t tG 


dr' d r >[r'g N (r')]H(r, r'). 


Note that the function H(r, r') can be shown to be proportional to the elliptic function K{x): 

H(r,r') = — K{ r —) 

7Tr > V > 


(17) 


(18) 


with r> (r<) the larger (smaller) of r and r'. 

The proof is quite straightforward. Since some properties of the Bessel functions are very important in the dynamics 
of the spiral disk, as well as many disk-like system, we will show briefly the proof for heuristical reason. One of the 
purpose of this derivation is to clarify that there are different definitions for the elliptic functions written in different 
textbooks. Confusion may arise applying the formulae in a wrong way. 

Note first that the Bessel function Jo(x) has an integral representation nm 

i r 

J 0 (x) = — cos[a:sin0]d0. 

n Jo 


( 19 ) 
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Given the delta function represented by the plane wave expansion: 

1 [°° 

S(x) = — / dkexp[ikx], 

27r J -oo 

One can show that, with the integral representation of Jo, 


dxcos[kx] Jo(x) = 


VT^k 2 


( 20 ) 


( 21 ) 


for k < 1. On the contrary, above integral vanishes for k > 1. Therefore, one can apply above equation to show that 


. 2 f */ 2 

dkJo(kx)Jo(k) = — 


dd 




x 2 sin 2 ( 


= —K{x) 

7r 


( 22 ) 


with the help of the Eq. m again. The last equality in above equation follows exactly from the definition of the 
elliptic function K. 

There is an important remark here. Note that the elliptic functions E(k ) and K(k) used in this paper, and Ref. 
j, 23] as well, are defined as 


r "V 2 , 

K(x) = / (l-x 2 sin 2 0)~ 1/2 de i (23) 

Jo 

p-K/2 

E(x) = / (l-x 2 sin 2 6») 1/2 J6» (24) 

J o 

which is different from certain textbooks. Some textbooks, and similarly some computer programs, prefer to define 
above integrals as K{x 2 ) and E{x 2 ) instead. In fact, one can check the differential equations satisfied by E(x) and 
K{x) that will be shown explicitly shortly in next section. Before adopting the equations presented in the texts, one 
should check the definition of these elliptic functions carefully for consistency. 

After a proper redefinition of k —> k'r' and write x = r/r' < 1, one can readily prove that the assertion m is 
correct. 

Note that the series expansion of the elliptic function K (x) is 


| E C ^ 2n (25) 

n=0 

with Co = 1 and C n = (2 n — l)!!/2n!! for all n > 1. Therefore, one can show that it agrees with the series expansion 

r 2 

of the hypergeometric function E(4, ±, 1; ^f), up to a factor tt/2. Indeed, one has: 


F(a, b 1 c; x 2 ) = ^ 

n —0 


( a)n(b)n 2 n 

(c)„n! 


(26) 


with {a) n = T(n + a)/T(a). It is straightforward to show that the expansion coefficients of Eq.s (1251) and (1261 agree 
term by term. Therefore, the derivation shown above agrees with the result shown in l23i l25l : 


H(r, r') = -F(l 1, 1; If) = —K(^). 

r> 2 2 r> n r> r> 


with r> (r<) the larger (smaller) of r and r'. 
As a result, one has 


/r(r) 


2t tG 


dAv 2 N {r')\ H(r, r')dr' 


(27) 


(28) 


given the identification g n = v %/ r ■ 
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III. SOME SIMPLE INTEGRABLE MODELS 


One can show that n( r ) is integrable with a Newtonian velocity described by 


v%(r) = 


y/r 2 + a 2 ’ 


(29) 


with Co and a some constants of parametrization. Note that the velocity function Ujv(r) vanishes at spatial infinity 
while r)jv(O) —> Cq with a non-vanishing value. One notes that the observed velocity at r = 0 is expected to be zero 
from symmetric considerations. We will come back to this point later for a resolution. Following Eq. o , one can 
write 


d( r ) = J kdkA(k)J 0 (kr) 

with A(fc) defined as 

m = r 

do 

To evaluate A (k), one will need the following formula: 


y/r 2 + a' 


- ; J\(kr)dr. 


dk exp[— kx\ Jo(fc) = 


VT 


(30) 


(31) 


(32) 


which follows from Eq. m by replacing x —> ix with the help of an analytic continuation. In addition, writing 
exp[— kx] dk as — d(exp[— kx\)/x and performing an integration-by-part, one can derive 


nOO 

/ dk exp[— kx] J\{k) = 1 — 

Jo 


y/1 + x 2 


(33) 


Here we have used the identity J' 0 = — Ji and the fact that Jo(0) = 1. With the help of above equation, one can show 
that 


nOO 

/ dk( 1 — exp[— 2k}) J\{kx) = 
Jo 


/ o xy/x 2 + 4 

Multiplying both sides of above equation with J\[k'x)xdx, one can integrate above equation and obtain 

2 1 — exp[— 2k] 


dxJi(kx) 


yjx 2 + 4 


(34) 


(35) 


with the help of the closure relation 0 . After a redefinition of parameters k —>■ ka/2 and x —> 2 r/a, one can show 
that 


drJi(kr) 


1 — exp[— ak\ 


y/r 2 + a 2 


and hence 


A (k) = (1 — exp[— ak]). 

rZ 

In additon, equation (El can also be written as 

r°° i 

/ dk exp[— kx] Jo(ka) = 


Jo 

after proper reparametrization. Therefore, one has 

~r2 


y/x 2 + a 2 


(36) 


(37) 


(38) 


1*00 

M(r) = J dk{l - exp[-ak])J 0 (kr) 


C 2 

°o 

2t tG 


y/r 2 + a 2 


( 39 ) 
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Hence this model is integrable as promised. Also, as mentioned earlier in this section, this model with a small Co and 
properly adjusted a can describe the velocity profile pretty nice in the case of MOND. 

One can eliminate the non-vanishing constant by two different methods. The first method is simply subtracting 
two integrable v n( t > C 0 , Oj), namely, define the new Newtonian velocity as 


v No( r ) — Q) 


a i 


a 2 


\Jr 2 + a\ yjr 2 + a 2 


(40) 


with Co and a\ > a 2 some constants of parameterizations. This new velocity function is also integrable due to the 
linear dependence of the function v %( r ) in Eq. G3- In addition, it vanishes at r = 0 and approaches 0 at spatial 
infinity r —> oo. Therefore, one has 




1 


-Jr 2 + d\ \Jr 2 + a\ 


(41) 


directly from Eq. In addition, one can also show that j2^| higher derivative models defined by 


V Nn{r) = -Cn(-i) 


_ l)!(2n 2 k)\ l- 2 fc/ 2 2 '-n+fc- 1/ 2 ( 2n )- / 2 . 2s-n+l/2 

+ ’ 2 2n n\ [ + ’ 

(42) 


da 2 ' Vr 2 + a 2 2 2n ~ 1 {k - l)!(n - fc)!(2Jfc - 1) 


are also integrable and give the mass density as 

'1 1 




Vr 2 


C 2 d 

w n /_ _ \r, 

27rG l da 2 ’ (Vr 2 +a 2 J ~~ 2nG 2 2n n\ 


G n ( 2n ) ! , 2 , a 2\-n- 1/2 

o2n„l v / 


(43) 


Mote that both v% n and fi n (f) are in fact functions of (r 2 + a 2 ) n 1 / 2 with appropriate combinations. For example, 
given the Newtonian velocity 


U N 1 


(r) = 


C 2 r 2 


a(r 2 + a 2 ) 3 / 2 
the corresponding mass density will be given by 


= C 2 


MiW = 


a(r 2 + a 2 ) 1 / 2 (r 2 + a 2 ) 3 / 2 

C 2 


2nG(r 2 + a 2 ) 3 / 2 


(44) 


(45) 


For convenience, we have absorbed a common factor 1/2 into C 2 . Note that vf^ 1 (r —* oo) —► 0 and v\(r = 0) = 0 in 
this model. 


A. Newtonian Model 


Consider the case of Newtonian model that observed velocity v and Newtonian velocity vjy are identical. This model 
is known to require the presence of dark matters [27). Since the velocity has to vanish at r = 0 and goes to a constant 
at spatial infinity, one will show that an additional constant term added to the vfj will provide both resolution at the 
same time. Indeed, another way to eliminate the non-vanishing velocity at r = 0 is to introduce a constant velocity 
by noting that 

1 f 00 C 2 

f ‘ (r) = 2^1 kdkA mM*r) = ^ (46) 

with A(fc) defined as 


Here we have used the identity 


A (k) = C 2 



J\{kr)dr 




dkJ n (kr ) 


1 


(47) 


r 


(48) 




































which follows directly from Eq. flSTCBI . 

Therefore, one can show that the Newtonian velocity given by the form 


7 


vl{r) = v% 0 (r) = Cl 1 


a 

yV 2 + a 2 


is induced by the mass density of the following form 


Mo (r) 


C 2 0 1 

2nG yjr 2 + a 2 


(49) 


(50) 


Here Co and a are some constants of parametrization. Note that the Newtonian velocity m vanishes at r = 0 and 
Vjv —► Co at spatial infinity as promised. Therefore, this velocity profile goes along with the observed RC of the 
spirals. Hence this model can be used to simulate the dynamics of spiral galaxies which requires the presence of dark 
matters. Hence we will call this model with v = vjyris the Newtonian model that normally requires the existence of 
dark matters. 

This zero mode is however the only known integrable modes for velocity profiles. Successive differentiating v 2 N in 
this model simply turns off the constant term Cg in Eq. Therefore, higher derivative modes derived from further 

differentiation of the Newtonian potential V N with respect to —a 2 , 


£ = ^(-3 a2 )"ug = C£(-3 oa ) r 


. vV 2 + a 2 . 


(51) 


will simply take away the leading constant term from the higher order modes. Therefore, velocity profiles v n (r) will 
be quite different from vq in this case. Note again that the velocity v n is induced by the corresponding mass density 


Mn(r) 


C 2 (2 n 1)!! , o , 2\ —n —1/2 
27rG 2™ [ + ’ 


(52) 


One of the advantage of these well-behaved smooth velocity functions is that they can be used as expansion basis 
for simulation of velocity profiles. Thanks to the linear dependence of the u^-(r) in the mass density function /z(r), 
one can freely combine any integrable modes of velocity to obtain all possible combinations of integrable models. For 
example, the model with 


v 2 {r) = v%(r) = v o( r , ° 0 i, Oj) = E c oi 
i,j i,j 



(53) 


is integrable and can be shown to be derived by the mass density 

MW = ^°( r ’ Co*, aj) = J2 o°r< r— — = ( 54 ) 

i i 2nG sj t 2 + a 2 

with C nl and cij all constants of parametrization. Here summation over i,j is understood to be summed over all 
possible modes with different spectrum described by Coi and a,. This velocity vanishes at r = 0 and goes to 

V 2 (r)-J 2 C oi ( 55 ) 

i,j 

at spatial infinity. One can also add higher derivative terms w 2 (r) to the velocity profile r> 2 (r). Since higher derivative 
velocity v n vanishes both at r = 0 and spatial infinity, these additional terms will not affect the asymptotic behavior 
of v 2 at r = 0. Therefore, one needs to keep at least one zeroth order term in order for v to be compatible with the 
RC data. 

To be more specific, one can consider the model 


\r) =vjf(r) = 


E 


vl{r,C n 




i,3 


1 - 


c 2 nk {-d a ,f 


k,l,n 


ai 


sjr 2 +i 


(56) 





























which be derived by the mass density 


K r ) = ^21 


i (^5 Cni ? dj ) — ^ ^ 


r 2 

L 0 i 


2ttG 


’ 2 + a 2 


E 


r 2 

^nk 


(2 n - 1)!! 


2ttG 


(r 2 + a 2 ) n 1/2 


(57) 


with C 0 i and Oj are all constants of parameterization. Here the summation over n is to be summed over all n > 1. 
The velocity of these models will vanish at r = 0 and goes to 

v 2 (r)^J2 C Oi (58) 

at spatial infinity. Therefore, these models turn out to be good expansion basis of any RC data for Newtonian 
dynamics. 

In practice, one may fit the v 2 in expansion of these modes in order to analyze the RC in basis of of these basis modes. 
This helps analytical understanding of the spiral galaxies more transparently. The properties of each modes is easy 
understand because they are integrable. The corresponding coefficients C n i and cij will determine the contributions of 
each modes to any galaxies. One will be able to construct tables for spiral galaxies with the corresponding coefficients 
of each modes. Hopefully, this expansion method originally developed in |23| will provide us a new way to look at the 
major dynamics of the spiral galaxies. 


B. Milgrom model 


For the case of Milgrom model (2), the Newtonian velocity vn and the observed velocity v are related by the 
following equation 


v 2 (r) = 


V^{r) + sjV^{r) + AV^{r)glr 2 


1/2 


Therefore, one can show that a galaxy with a rotation curve given by 


C A n 2 

4/ \ '-'0 U 

W) = 


1 + \ 1 + 


4(/gr 2 (r 2 + a 2 ) 


2 (r 2 + a 2 ) 

is the corresponding velocity induced by the Newtonian velocity 

v 2 m ( r ) = 


cy 


C 2 a 


(59) 


(60) 


(61) 


\Jr 2 + a 2 

Here Co and a are some constants of parametrization. Therefore, the corresponding mass density is given by Eq. (® 

, x Co r i i 
Mr) = 


2t tG 


r vV 2 + a 2 


Note that Vq approaches a constant goC^a at spatial infinity. This asymptotically flat pattern of the velocity is 
compatible with many observations of the spirals. This Newtonian velocity does not, however, vanish at the origin. 
Indeed, one can show that u(0) —> Go 7 ^ 0 Therefore this would not be a good expansion basis for the physical spirals. 
The non-vanishing behavior of v at r = 0 can be secured by considering the refined model (ECT : 


v No( r ) — Cq 


a 1 


a 2 


y \Jr 2 + a\ y/r 2 + a\ J 

with Gq, ai > 02 some constants of parameterizations. The corresponding velocity function v can be shown to be 


si 4 

4/ \ °0 

v o( r ) = -7T 



2 

f 


1 / 2 ) 

fll 02 


1+ 

X| 4 5 2 r 2 (r 2 +a 2 )(r 2 +o 2 ) 


\/r 2 + a\ yjr 2 +a\ 

1 


Cq (aiyV 2 a 2 y/r 2 + afj 

1 


(62) 
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This new velocity function v is hence induced by the Newtonian velocity cm . In addition, Note that v 4 approaches 
a constant <7o^o( a i — ^ 2 ) spatial infinity and vanishes at r = 0. This agrees with the main feature of the observed 
asymptotically flat rotation curve of many spirals. As a result, the corresponding mass density of this model is given 
by Eq. 14II) : 


w(r) = 


\fr 2 + a 2 \fr 2 + a\ 


In addition, a model with a velocity, in the case of Milgrom model, of the form 


v, = 


Cfr 4 + \J Cfr 8 + ACfi' 6 gQa 2 (r 2 + a 2 ) 3 


2 a 2 (r 2 4- a 2 ) 3 

can be shown to be induced by the Newtonian velocity of the following form given by Eq. Cl : 

r 2 r 2 

2 / \ ( ~'l r 

V N lW = , 2 . 2 N 3/2 - 

ayr 2 + a 2 ) 6 ' 2 

Therefore, this model is derived by the mass density C51) : 


(63) 


C 2 

= 2 7rG(r 2 + a 2 ) 3 / 2 ' 

Note that v 2 (r —> oo) —> C\\]go/a and v\(r = 0) = 0 in this Milgrom model. In addition, the corresponding 
Newtonian model also has the same properties: Ujvi(r) vanishes both at r = 0 and r —* oo. Therefore, this model 
also appears to be a realistic model in agreement with the asymptotically flat rotation curve being observed. 

Note again that further differentiation of the Newtonian velocity v% with respect to —a 2 will derive integrable 
higher derivative models: 


v Nn 


(%(-d a ,y 


a 

\Jr 2 + a 2 


(64) 


Therefore, this velocity will be derived by the mass density 


Vn{r) 


Gn ( 2n 1 )‘‘ /y2 , 2\ —n —1/2 

27 tG 2™ [ + ’ 


(65) 


One of the advantage of these well-behaved smooth velocity functions is that they can be used as expansion basis 
for simulation of velocity profiles. Thanks to the linear dependence of the uXr(r) in the mass density function /z(r), 
one can freely combine any integrable modes of velocity to obtain all possible combinations of integrable models. For 
example, the model with 


U N 


(r) = ^no( r, C 0 i, aj,bj) = Y C < 


\J r2 + a2 j \/ r2+b l 


is integrable and can be shown to be derived by the mass density 


M r ) = = 5Z 


r 

'-'I 


0 i 


27 tG 


v' r2 + i> “ 


( 66 ) 


(67) 


with C n i and aj all constants of parametrization. The velocity also vanishes at r = 0 and goes to 

vl(r) (68) 

at spatial infinity. This will in turn make the corresponding observed Milgrom velocity v 2 approaches the asymptotic 
velocity v ^ j C^go^aj — bj )] 1 ^ 2 . One can also adds higher derivative terms v% n (r) to the velocity profile v%{r). 
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Since higher derivative velocity v n goes to zero faster than the zero-th derivative term at spatial infinity, these adding 
will not affect the asymptotic behavior of v 2 at spatial infinity. Therefore, leading order terms will determine the 
asymptotic value of v. To be more specific, one can consider the model 

v %( r ) = Cni ' - 5Z' 

hj 

which be derived by the mass density 

^0) = r ’ a j) = g 

i i,j 

with Coi and a,j all constants of parametrization. Note again that the velocity v 2 of these models also vanishes at 
r = 0 and goes to 

(71) 

i,3 


' 2 + a2 j 


b l 


+ '£c 2 nk (-d a *y 


k,l‘ 


ai 


sjr 2 + af 


(69) 


r 2 + b 2 


r 2 + a 2 


V" (2n 1 ) ! '( r 2 +a 2 )—-1/2 


22^ 27 tG 2” 

n,k,l 


(70) 


corresponding to 


2 (r) 


go 


Cq i( a j bj) 


7J 


1/2 


(72) 


at spatial infinity. Therefore, these models turn out to be good expansion basis for v%(r) of any RC data for Milgrom 
models. 

In practice, one may convert the RC data from v to vn following Eq. m and then fit the resulting v ^ in expansion 
of these modes in order to analyze the RC in basis of of these basis modes. This helps analytical understanding of the 
spiral galaxies more transparently. The properties of each modes is easy understand because they are integrable. The 
corresponding coefficients C n i and a,j will determine the contributions of each modes to any galaxies. One will be able 
to construct tables for spiral galaxies with the corresponding coefficients of each modes. Hopefully, this expansion 
method originally developed in 12/1 will provide us a new way to look at the major dynamics of the spiral galaxies. 


C. Famaey and Binney model 

For the case of FB model ©, the Newtonian velocity vn and the observed velocity v are related by the following 
equation 


v 2 = 


i/ vn 4, + 4 gorvjf + v 


(73) 


Therefore, one can show that a galaxy with a rotation curve given by 


2 / , C 2 a 
v 0 {r) = 


2 sjr 1 + a 2 


1 + 


1 + 


4:g 0 rV r 2 + a 2 


1 / 2 ' 


(74) 


is the corresponding velocity induced by the Newtonian velocity 


v 2 m{r) = 


\Jr 2 + a 2 


(75) 


Here Cq and a some constants of parametrization. Therefore, the corresponding mass density is given by Eq. 01 


r 2 

Mr) = ° 


1 


1 


2itG \_r yj r 2 + a 2 . 
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Note that Vq approaches a constant goC^a at spatial infinity. This asymptotic flat pattern of the velocity is compatible 
with many observations of the spirals. This Newtonian velocity does not, however, vanish at the origin. Indeed, one 
can show that uo(0) —> Co 4 1 0- Therefore this would not be a good expansion basis for most physical spirals. The 
non-vanishing behavior of v at r = 0 can be secured by considering the refined model BUI) : 


v No( r ) — Q) 


a i 


a 2 


\Jr 2 + a\ \Jr 2 + i 


with Cq, a\ > a 2 some constants of parametrization. The corresponding velocity function v can be shown to be 


r 2 

L o 

<2l 

a 2 

f + 

x j 4 ff0 r(r 2 +a 2 1 ) 1/2 {r 2 + cl%) 1/2 

1/2' 

2 

\Jr 2 + a\ 

yV 2 + a\ 

Cq (yii\/ f 2 +ai~ a 2 yjr 2 + afj 

> 


(76) 


This new velocity function vo is hence induced by the Newtonian velocity BUI . In addition, Note that Vq approaches 
a constant goC^ai — a 2 ) spatial infinity and vanishes at r = 0. This fits the main feature of the asymptotically 
flat rotation curve of the spirals. As a result, the corresponding mass density of this model is given by Eq. Gl: 


'‘" W = Sc 


1 


1 




In addition, a model with a velocity, in the case of FB model, of the form 

Cfr 2 [ 4 g 0 a(r 2 + a 2 ) 3 / 211 / 2 


2 a(r 


2j3/2 


1 + 


1 + 


C 2 r 


can be shown to be induced by the Newtonian velocity of the following form given by Eq. BUI : 

C 2 r 2 

2 ( \ u l r 

^Nlv) / 2 i 2'\3/2 ’ 

a(r z + a z ) 6 ' z 

Therefore, this model is derived by the mass density 63): 


(77) 


r 2 

( \ __w_ 

^ ll[r> - 2tt G{r 2 +a 2 ) 3 / 2 ' 

Note that vf(r —> oo) —> C\\Jgo/a and v\{r = 0) = 0 in this FB model. The corresponding Newtonian model also 
has the same limit: vni{t) goes to 0 in both r = 0 and r —> oo limits. Therefore, this model appears to be a more 
realistic model compatible with the flat rotation curve being observed. 

Note again that further differentiation of the Newtonian velocity vjj with respect to — a 2 will derive integrable 
higher derivative models: 


v 


Nn 


= c 2 n (-d a ,y 


Vr 2 + 


a* 


Therefore, this velocity can be shown to be derived from the mass density 

^2 


Mr) = 1 )!! (r 2 +o 2 )-"- 1/2 - 


2t tG 


(78) 


(79) 


One of the advantage of these well-behaved smooth velocity functions is that they can be used as expansion basis 
for simulation of velocity profiles. Thanks to the linear dependence of the v n( v ) in the mass density function g(r), 
one can freely combine any integrable modes of velocity to obtain all possible combinations of integrable models. For 
example, the model with 


v n ( r ) — ( r > Co i i Uj ,bj) = Cq^ 


ij 






(80) 
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is integrable and can be shown to be derived by the mass density 


MO) = E «>( r > C °” °J') - E 


r‘ 


0 i 


27 tG 


y / r2+b2 


3 J 


(81) 


with C n i and aj > bj are all constants of parameterizations. The velocity will vanish at r = 0 and goes to 


(82) 


h3 


at spatial infinity. This will in turn make the corresponding observed FB velocity v approach the asymptotic velocity 
v oo Ej j Coi5o(ai — bj)] 1 / 2 . One can also adds higher derivative terms V Nn( r ) to the velocity profile Since 

higher derivative velocity v n goes to zero faster than the zero-th derivative term at spatial infinity, these adding will 
not affect the asymptotic behavior of v 2 at spatial infinity. Therefore, the leading order terms will determine the 
asymptotic behavior of the RC. 

To be more specific, one can consider the model 


v 2 N (r) = J2 v U r , C ni , Oj) = ^2C{ 
which be derived by the mass density 


h3 


r 2 + a 2 Jr 2 + b) 


E CnkhdafY 

n,k,l 


ai 


V r2 + a i 


(83) 


/ i ( r ) = E ^”( r ’ c ni,aj) = E 


r 


0 i 


2t tG 




E S ,2 ”.„ 1)!! (-- + (84) 


„k,l 


2ttG 2 11 


with Coi and aj all constants of parameterization. The velocity of these models vanishes at r = 0 and goes to 


U N 


H-E^ b ‘ 


(85) 


h3 


corresponding to 


- 1/2 

90^2 C 0i( a 3-hj) ( 86 ) 

i,j 

at spatial infinity. Therefore, these models turn out to be good expansion basis for v\, of any RC data for FB models. 

In practice, one may convert the RC data from v to vn following Eq. (17.31) and then fit the resulting vj^ in expansion 
of these modes in order to analyze the RC in basis of of these basis modes. This helps analytical understanding of the 
spiral galaxies more transparently. The properties of each modes is easy understand because they are integrable. The 
corresponding coefficients C n i and aj will determine the contributions of each modes to any galaxies. One will be able 
to construct tables for spiral galaxies with the corresponding coefficients of each modes. Hopefully, this expansion 
method originally developed in t23l will provide us a new way to look at the major dynamics of the spiral galaxies. 



IV. COMPACT AND REGULAR EXPRESSION 

In order to put the integral in a numerically accessible form, Eq. for /i(r) can be written as a more compact 
form with a compact integral domain x £ [0,1]: 


mM = 


1 


7 r 2 Gr 

1 r 1 r 

d x [v 2 N (rx)]K(x)dx- / d v [v%(-)\ K(y)ydy 

uo J o y 


(87) 
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Here one has replaced x = r'/r and y = r/r' in above integral. One of the advantage of this expression is the numerical 
analysis involves only a compact integral domain 0 < x < 1 instead of a open and infinite domain 0 < r —> oo domain. 
Even most integral vanish quickly enough without bothering the large r domain, the compact expression will make 
both the numerical and analytical implication more transparent to access. 

Note that the function K(x) diverges at x = 1. It is, however, easy to show that K(x)dx —■> 0 near x = 1. Usually, 
one can manually delete the negligible integration involving the elliptical function K(x —> 1) to avoid computer 
break-down due to the apparent singularity. 

It will be, however, easier for us to perform analytic and/or numerical analysis with an equation that is free of any 
apparently singular functions in the integrand. Indeed, this can be done by transforming the singular elliptic function 
K{x) to regular elliptic function E{x). The advantage of this transformation will be used to evaluate approximate 
result in the next sections. Therefore, one will try to convert the apparently singular function K{x) into a singular 
free function E{x) by performing some proper integration-by-part. 

One will need a few identities satisfied by the elliptic functions E and K. Indeed, it is straightforward to show 
that E(x) and K(x) satisfy the following equations that will be useful in converting the integrals into more accessible 
form: 


x{l — x 2 )K"(x) + (1 — 3x 2 )K'(x) — xK(x) 

= 0, 

(88) 

x{l — x 2 )E"{x) + (1 — x 2 )E\x) + xE(x) 

= o, 

(89) 

xE'(x) + K{x) = E{x ), 


(90) 

E\x) + ~^E{x) = K’(x), 

1 — x z 


(91) 

xK'(x) + K(x) = -— ^—=E(x). 

1 — x z 


(92) 


Note that Eq. (HJUll can be derived directly from differentiating the definition of the elliptic integrals Eq.s (1231241) . In 
addition, Eq.s (1911921) can also be derived with the help of the Eq. 189f) . 

When one is given a set of data as a numerical function of Ujv(rx), it is much easier to compute dv^^k = rx)/dk 
instead of d x vjv(rx). Therefore, one will need the following converting formulae: 

d x [vjf(rx)] = r[v%(rx)]', 


d v [v 2 N (-)] = -^[v 2 N (-))'. 


Therefore, one is able to write the Eq. (IH7I) as: 


K r ) = 


r 2 G 


dv 2 N {rx) K{x)dx 


/ dv 2 N (-)I<(y)y i dy 
/ o V 


(93) 


Here dv^(r) = d r Kr(0] = [<M]' with ' denoting the differentiating with the argument r, or rx. 

Note that the part involving the integral with dvj^(r/y) is related to the information in the region r' > r. Here r 
is the point of the derived information such as y(r), and r' is the source point of observation v(r') in the integrand. 
Therefore, this part with source function r/y contains information exterior to the target point r. On the contrary, 
the source term with function of rx represents the information interior to the target point r. 

Most of the time, the source information beyond certain observation limit r' = R becomes un-reliable or unavailable 
due to the precision limit of the observation instruments. One will therefore need to manually input the missing data 
in order to make the integration result free of any singular contributions due to the boundary effect. One will come 
back to this point in section IV. 

In addition, Eq. (PI) can be used to derive the total mass distribution M (r) of the spiral galaxy via the following 
equation: 


M(r) = 2n f r'dr 1 (94) 

Jo 

Note that the velocity function vn shown previously in this paper is the rotation velocity needed to work with the 
total mass of the system in the Newtonian dynamics. One can derive the velocity v(r) that works with the dynamics 
of MOND with the relation given by Eq. © and @. In short, v —> vn in the limit of the Newtonian dynamics. 
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Throughout the rest of this paper, we will discuss the application of these formulae both in the case of the Newtonian 
dynamics with dark matter and in the case of MOND. Therefore, we will first derive the velocity function VRr(r) from 
v(r) in the case of MOND. As mentioned above, they follows the relation given by Eq. m and 0. 

Two different models will be studied later: 

Case I: Milgrom model 

Indeed, if one has g(r) = v 2 (r)/r in the case of Milgrom model, one can show that 


v 2 N (r) 


v 4 (r) 

\Jv 4 {r) +g%r 2 


(95) 


In dealing with the exterior part involving ro/y, one has to compute dvj^(r) at large r. By assuming v(r) —> vr = 
v(r = R), one can show that: 


dv 2 N (r > i?) —> — - 


4 2 

Vr 9 o r 


v r + 9o r2 ) 3/2 ' 


(96) 


Here R is the radius of the luminous galactic boundary. Mostly, the flatten region of RC becomes manifest beyond 
r > R. In addition, vn(tx) and VRi{r/y) take the following form: 


v 2 n (rx) 


v 4 {rx) 

\/v 4 {rx) + gQr 2 x 2 ’ 


(97) 


<0 = 


l (r/y)y 


9 \fv 4 {r/y)y 2 + g%r 2 

Therefore, the surface mass density y,(r) can be written as, with the velocity vr r(r) given above, 

1 


mW = 


2t tG 


dr'\ 


\Jv 4 T g^r ' 2 


H{r , r')dr' 


(98) 


(99) 


in the case of Milgrom model. 

Case II: FB model 

Similarly, if one has g(r) = v 2 (r)/r in the case of FB model, one can show that 


v%(r) = 


v 4 (r) 


v 2 (r) + g 0 r 

By assuming v(r) —> vr = v(r = R), one can also show that: 


dv%(r > R) —> — 


vfigo 


(vr + 9or) 2 


In addition, ujv(rx) and visr(r/y) take the following form: 


v 2 N {rx) = 


v 4 (rx) 


v 2 (rx) + gorx ’ 


( 100 ) 


( 101 ) 


( 102 ) 


v 2 n (~) = 


\r/y)y 


V v 2 (r/y)y + g 0 r' 

Therefore, the surface mass density y,(r ) can be written as, with the velocity u/v(r) given above, 


= if d A-^TZ - 7 ] H (r,r')dr’ 

2t tG J o v 2 + g 0 r' 


(103) 


( 104 ) 


in the case of FB model. We will try to estimate the exterior contribution of these two models shortly. 
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As mentioned above, it is easier to handle the numerical evaluation involving regular function E(x) instead of the 
singular function K{x). Therefore, one can perform an integration-by-part and convert the integral in Eq. 18711 into 
an integral free of singular function K(x). The result reads, with fi(r) = @(r)/ (n 2 Gr), 


e(r)= f E'(x)Vi(rx)dx 

Jo 1 x Jo 


/ 0 x A Jo 

The last term on the right hand side of above equation can be integrated by part again to give 

r± V$(r/ir) ~ xV${rx) 


0(r) = 


/o 


1 — X 2 


-E(x)d: 


x 


+ [ E{x)d x V^{rx)dx - V$(r). 

Jo 


Hence one has 


H(r) = 


7r 2 Gr J o 


Vn(t/x) - xV£(rx) 
l~x 2 


E{x)dx 


E(x)d x Vj$(rx)dx - V£{r)]. 


(105) 


(106) 


(107) 


Note again that the integral involving Vn(t'x) carries the information r' < r while the integral with vn(j/x) represents 
the contribution from r' > r by the fact that 0 < x < 1. As promised, one has transformed the singular I\ function 
into the regular E function. 

Although there are still singular contribution like 1/(1 — x 2 ) in the integrand, it is easier to handled since we know 
these functions better than K function. This is because we only know a formal definition of this function via a set 
of definitions. Even we have a rough picture about the form of K(x) and K'{x). Numerical and analytical analysis 
could be difficult as compared to dealing with the more well-known function like 1/(1 — x 2 ). We will show explicitly 
one of the advantage of this equation in next section when one is trying to estimate the contribution from a model 
describing the missing part of observation. 


V. CONTRIBUTION FROM THE ASYMPTOTIC REGION 


In practice, measurement in far out region is normally difficult and unable to provide us with reliable information 
beyond the sensitivity limit of the observation instrument. One often can only measure energy flux and rotation curve 
within a few hundred kpc from the center of the galaxy. Beyond that scale of range, signal is normally too weak to 
obtain any reliable data. Therefore, one has to rely on various models to interpolate the required information further 
out. 

It it known that, contrary to the spherically system, exterior mass does contribute to the inner region. Therefore, 
it is important to estimate the exterior contribution carefully with various models. In this section, we will study a 
velocity model with a flat asymptotic form and its contribution to the inner part in both Newtonian dynamics and 
MOND cases. One of the purpose of doing this is to demonstrate the advantage of the regular function formulae one 
derived earlier. 

For a highly flatten galaxy, formulae obtained earlier in previous section has been shown to be a very useful tool 
to predict the dynamics of spiral galaxies. It is also a good tool for error estimation. Possible deviation due to the 
interpolating data can be estimated analytically more easily with the equations involving only regular elliptic function 
E(x). Note again that another advantage of Eq. (111)71) is that the r-dependence of the mass density n{r) has been 
extracted to the function vn- This will make the analytical analysis easier too. 

Assuming that the observation data v is only known for the region r < R 1 the following part of Eq. (117771) 


5p(r < R) = 


1 

n 2 Gr 



Vn^/x) 

1-x 2 


E(x) 


(108) 


represents the contribution of p(r < R) from the unavailable data VN(r') beyond the point r = R. To be more 
precise, one can write n(r < R) = + 5p(r) with the mass density /r<(r) being contributed solely from the 
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vj v(r') = V]y(r/x) data between 0 < r' < R or equivalently r/R < x < 1. Explicitly, n<(r) can be expressed as 





(109) 


As mentioned above, one will need a model for the unavailable data to estimate the contribution shown in Eq. 
TO. We will show that a simple cutoff with v(r > R) = 0 will introduce a logarithmical divergence to the surface 
density /i(r = R). The divergence is derived from the singular denominator 1 — x 2 in Eq. Ill (171) . The factor 1/(1 — x 2 ) 
diverges at x = 1 or equivalently r' = r. A smooth velocity function Vn{t) connecting the region R — e < r < R + e is 
required to make the combination [V^(r/x) — xV// r (rx)\/( 1 — or) regular at x = 1. Here e is an infinitesimal constant. 

Case I: Milgrom model 

Let us study first the case of Milgrom model. Since most spiral galaxies has a flat rotation curve v(r R) —> vr 
with a constant velocity vr. For our purpose, let us assume that vr = v{r = R) for simplicity. Hence the deviation 
(fTTTHl) can be evaluated accordingly. Note again that this simple model agrees very well with many known spirals. 

Note first that the Newtonian velocity v m (r) is given by Eq.s and (HQ) with v(r > R) = Vr. After some 
algebra, one can show that 


Shm(t < R) = 


Em 


r/R 


r v 2 N (r/x) 

1 - X 2 


n 2 Gr 


dx 


t 2 Gr 


rr/R 

/ dx 

'v%{r/x) ' 

1 0 

1 — X 2 


E m Vr 


■K 2 Gr y /vj t + g%r 2 


In 


r V v R + 9 qE 2 + R-Vvr +9or 2 

(, gor + V v r + 9%r 2 )VR 2 - r 2 


E (x) 


( 110 ) 

( 111 ) 


( 112 ) 


Note that tt/ 2 > E(x) > 1 is a monotonically decreasing function with a rather smooth slope. The rest of the 
integrand is also positive definite. Therefore, one can evaluate the integral by applying the mean value theorem for 
the integral ifTTUl) with Em = E(x = Xm) the averaged value of E(x) evaluated at xm somewhere in the range 
0 < xm < 1 . 

Case II: FB model. 

Let us study instead the case of FB model. Let us also assume that vr = v(r = R) for simplicity. Note first that 
the Newtonian velocity vn(v) is given instead by Eq.s cm and cm with v(r > R) = Vr. After some algebra, one 
can show that 


with 


< R) = 


Er 

TT 2 Gr 



v 2 N( r /x) 

1 — X 2 


efv 4 r r 
2 Gr J 0 


1 

.(! - x 2 )(v 2 r + gor /x) _ 


Efv% t 

n 2 Gr 


(113) 


I = 



X 

(' l-x 2 )(x + g 0 r/v 2 R ) _ 


(114) 


Note that n/2 > E(x) > 1 is a monotonically decreasing function with a rather smooth slope. The rest of the integrand 
is also positive definite. Therefore, one can evaluate the integral by applying the mean value theorem for the integral 
cm with Ep = E{x = xp) the averaged value of E(x) evaluated at Xf somewhere in the range 0 < Xf < 1. After 
some algebra, one can show that 


S/j,F(r < R) = 


EfVr 

2iT' 2 Gr 


9 -ln(l + r/R) --ln(l - r/R) - 

v 2 R -g 0 r v R + gor 


2gor 


4 2 9 

v r - 9o r- 


ln[l + v%/(g 0 R)\ 


(115) 


Case III: Newtonian case. 
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Similarly, one can also evaluate the mass density in the case of Newtonian dynamics with dark matter in need. Let 
us assume Vn{t > R) = Vr for simplicity again. As a result the deviation of mass density required to generate the 
rotation curve vn(t) = v(r) can be shown to be 


8g n{ r < R) 


v 2 N (r/x) 
1 — x 2 

e nv\ ln R + r 
2n 2 Gr R — r 


En 

n 2 Gr 


r r /re 


dx 


(116) 


after some algebra. Note that, in deriving above equation, one also applies the mean value theorem with En = E(x = 
xjv) the averaged value of E(x) evaluated at Xn somewhere in the range 0 < Xn < 1. 

To summarize, one has shown clearly with a simple model for the asymptotic rotation velocity that formulae with 
regular E function appears to be easier for analysis. This is mainly due to the fact that E{x) is smoothly and 
monotonically decreasing function within the whole domain x £ [0,1]. In most cases, mean value theorem is very 
useful in both numerical and analytical evaluations. 

In addition, one notes that the mild logarithm divergent terms appeared in the above final results are due to the 
cut-off at r = R. A negative and equal contribution from the interior data will cancel this singularity at r = R. 
To be more precisely, if we turn off the v function abruptly starting the point r = R by ignoring the exterior region 
contribution, a logarithmic divergence will show up at r = R accordingly. The appearance of the logarithm divergence 
also emphasize that the boundary condition of these physical observables at r = R should be taken care of carefully 
to avoid these unphysical divergences. In practice, one normally adds a quickly decreasing Vn{x > R) to account 
for the missing pieces of information and to avoid this singularity. Numerical computation may, however, bring up a 
small peak near the boundary r = R if the matching of Vn at the cutoff is not smooth enough. One will also show 
that similar singular behavior also appears in the final expression of the velocity function derived from a given data of 
mass distribution in next section. Evidence also shows that exterior contribution should be treated carefully in order 
to provide a meaningful fitting result. 

In order to compare the difference of Sg for different models, we find it is convenient to write B = v^/^goR) and 
s = r/R such that A and r' both become dimensionless parameters. Note that B ~ 1.1 if we take vr ~ 250 km/s 
and R ~ 4.97 x 10 4 ly from the data of Milky Way. Therefore, B ~ 1.1 is typically a number slightly larger than 1. 
Hence, one can write above equations as 


Sy M (s < 1) 

5hf{s < 1 ) 
Sfj, N {s < 1) 


EMB 2 g 0 ^ 
t t 2 Gs^B 2 + s 2 
E F B 2 g 0 Tn(l + s ) 
2n 2 Gs B — s 


sV B 2 + 1 + V B 2 + s 2 
(s + \JB 2 + s 2 )V 1 — s 2 
ln(l — s) | 2sln(l + il) 

B + s + s 2 — B 2 ’ 


E N Bg 0 1 + s 
2tt 2 Gs 1 — s 


(117) 

(118) 
(119) 


In addition, one can estimate the deviation Sg at small s where s <C 1, or r <C R. The leading terms read: 


£>/iM(s < 1) 

8/j,f(s 1 ) 

SgN{s « 1) + O(s). 


Em go 
7 t 2 G 
Ef9o 
t t 2 G 


\/B 2 + 1 


1 + O(s), 


[B — \n(B + 1)] + O(s), 


ENgo 
7 t 2 G 


B 


( 120 ) 

( 121 ) 

( 122 ) 


at small r. Note that go/G ~ 0.18. Therefore, the most important contribution from the exterior contribution is 
near the boundary at r = R. Our result shows that special care must be taken near the boundary of available 
data. Appropriate matching data beyond this boundary is needed to eliminated the naive logarithm divergence. The 
compact expression also made reliable estimation of the deviation possible. 


VI. GRAVITATIONAL FIELD DERIVED FROM A GIVEN MASS DENSITY 

One can measure the flux from a spiral galaxy and try to obtain the mass density with the M/L = constant law 
f24ll . Even the M/L law is more or less an empirical law, it does help us with a fair estimate of the mass density 
distribution. We will focus again on the physics of a highly flattened spiral galaxy. Once the mass density function 
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is known, for the range 0 < r < R, one can also compute the gravitational field 3 j\r(r) from a given mass density 
function g(r). Once the function gN{ r ) is derived, one can derive g(r) following the relation given by Eq.s 0 and 
0 in the case of MOND. 

Indeed, one can show that < 7 Ar(r) is given by 


/•oo /»oo 

3 ,/v(r) = 27rG / kdk / r'dr' n(r')Jo(kr')Ji(kr) 

Jo Jo 

from the Eq.s 0 and O- Therefore one has 

/•OO 

g N {r) = 2TvG r’dr' n(r')H 1 (r, r') 

Jo 


with 


/•oo 

H\(r, r') = / kdkJ±(kr) J 0 (kr') = — d r H(r, r'). 
Jo 


By differentiating Eq. m, one can further show that Hi[r,r') becomes 


Hi(r, r') = 


- 




for r < r' 


n (r 2 — r' 2 ) y r 


E (—) for r > r '. 


(123) 

(124) 

(125) 

(126) 
(127) 


Note that one has used the differential equations obeyed by E and K shown in Eq.s Following the method 

shown in section IV, one can rewrite the equation as 


( r ) =4 G f dx 
Jo 


. ,xE(x) . , ,K(x ) 

^( ra )l^2 +^ r / x )\-^r 


E{ a 


c 2 (l — X 2 ) 


(128) 


after some algebra. In order to eliminate the singular function K ( x ), we can also convert K ( x ) into regular function 
E(x) following similar method. The result is 


gN(r) = 4G f dx 
Jo 


(129) 


g,(rx) 


xE(x) 


- H{r/x)[ 


E(x) E'{x ) 


+ 


1 — x 2 y 1 — X 2 

With an integration-by-part, one can convert E'[x) to a regular function E{x). The result is 


g N {r) = 2G/i*(r) — 27 t Gg,(r) + 4G f dx 

Jo 


(130) 


E{x) 


yr 2 (l — x 2 ) 


[x 3 fi{rx) — fi(r/x)\ — 


E[x)rg! (r/x) 


Here /r*(r) = lirn^^o 2 g{r/x)/x. If /.i(r —> oo) goes to 0 faster than the divergent rate of r, /r*(r) will vanish or remain 
finite. For example, if g{r > R) = 2 fj,RRr/(R 2 + r 2 ), one can show that /r*(r) = 4 unR/r. Here gn = /x(r = R). We 
will be back with this model in a moment. 

Similar to the argument shown in section IV, one can show that the terms with g,(r/x) in Eq. 111301 1 will contribute 


i-r/R 

Sg N o(r < R) = —4G / dxE(x) 
Jo 


[i[r/x ) rg!(rjx) 


\l-x 2 ) 


(131) 


to the function 3 Ar(r) due to the unavailable data g(r > R). Note, however, that part of the exterior region contribution 
has been evaluated via the integration-by-part process in deriving Eq. 0(1 involving vn{t/x ) andlU^ir). Therefore, 
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one should start with the complete version Ml in evaluating the deviation of giy due to the exterior part. Hence 
one should have 


nr / R 

Sg N (r < R) = —4 G / dxg(-) 
J n x 


E(x) 
1 — x 2 


E'{x) - 

x 


(132) 


For simplicity, we will assume the following form of mass distribution in the region r > R, 

. . 2 Rr 

g{r > R)= HR R2 + r2 (133) 

with Hr = g(r = R). Note that continuity of g(r = R) across the matching point r = R is managed to remain valid 
in this model. Moreover, g(r —> oo) —> 0 is expected to hold for the luminous mass of the spiral galaxies. After some 
algebra, one can show that 


SgN(r) 


R 

—4 Gg,R — 
r 


~ n + A ( r ) 


(134) 


with A(r) given by the integral 
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2b 


_( 1 ~y)( 1 + by) (l + by) 2 _ 

, 1 2b 

ay 


e(Vv) 


(l-y){l + by) (1 + by) 2 _ 


(135) 

(136) 


with the first two terms in Eq. (THHI the contribution from integration-by-part of E'{x). Here one has defined y = x 2 
and b = R 2 /r 2 for convenience. 

In addition, one also applies the mean value theorem to the integral 113511 by noting again that (i) 7r/2 > E(x) > 1 
is a monotonically decreasing function with a rather smooth slope, and (ii) the integrand is a positive function 
throughout the integration range. Therefore, the integral A(r) can be evaluated by applying the mean value theorem 
with Ei = E{x = X\) the averaged value of E(x) evaluated at x\ somewhere in the range 0 < Xi < 1. 

The remaining integral in A can be evaluated in a straightforward way and one finally has 


SgN{r) = 

. „ R 

-^Gg,R— 


Ei[l 


R 2 




(137) 


This is the contribution to the Newtonian field giv(r < R) due to the unknown exterior mass contribution. Note 
that similar singularity also appears at Sgj\r(r = R) which means that a careful treatment in modelling the unknown 
exterior mass is needed. 

Case I: Milgrom model. 

For the case of Milgrom model, one can show that 


g(r) = 


'glr{r) + \/g%(r) + 4:g 2 N (r)g { 


Therefore, the deviation 5g 2 (r) can be shown to be 


Sg 2 (r) = 


9 2 N {r) + 2<?o 


2v / ffAf( r ) + 4 fl , w( r )ffi 


Sg%(r) 


by solving the algebraic equations 0 - 0 . 
Case II: FB model. 

For the FB model, one can show that 


(138) 


(139) 


\J gN 2 + ^gogN + gN 

g = - 
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( 140 ) 
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Therefore, the deviation 5g 2 (r ) can be shown to be 


Sg(r) 


1 

2 


gjy(r) + 2g 0 

\/3w( r ) + 4gw(r)g 0 




(141) 


by solving the algebraic equations dust . 

Therefore, the deviation 5g{r) can be evaluated from above equation to the first order in <5gvr(r) with all gAr(r) 
replaced by gN<{f )■ Here gN<{r) is defined as the contribution of interior mass to the Newtonian field gw(r), namely, 


g N {r <R) = gN<{r ) + Sg N {r). 


(142) 


To summarize again, one has shown that the formulae with regular E function appears to be helpful in deriving 
the gravitational field strength for numerical and analytical purpose. 


VII. CONCLUSION 

We have reviewed briefly how to obtain the surface mass density fj,(r) from a given Newtonian gravitational field 
gN with the help of the elliptic function K(r) in this paper. The integral involving the Bessel functions is derived in 
detailed for heuristical reasons in this paper too. In addition, a series of integrable model in the case of MOND is 
also presented in this paper. 

One has also converted these formulae into a simpler compact integral making numerical integration more accessible 
and analytical estimate possible in this paper. The apparently singular elliptic function K(r) is also converted to 
combinations of regular elliptic function E(r) by properly managed integration-by-part. 

As a physical application, one derives the interior mass contribution /i(r < R) from the possibly unreliable data 
v(r > R) both in the cases of MOND and in the Newtonian dynamics. Detailed results are presented for Milgrom 
model and FB model for the case of MOND. In particular, the singularity embedded in these formulae are shown to be 
a delicate problem requiring great precaution. Similarly, one also tries to derive similar results for the corresponding 
gN from a given /r(r). One also presents a simple model of exterior mass density /i(r > R) as a simple demonstration. 
The corresponding result in the theory of MOND is also presented in this paper. 

In section III, one has studied many solvable models in details. Analysis is generalized to Newtonian models, 
Milgrom models as well as the FB models. In practice, may convert the RC data from v to Vn following the 
transformation formula of either Milgrom model or FB model. For the case of Newtonian model, the RC data 
gives exactly the Newtonian velocity, namely, v = ujv- One can then fit the resulting v % in expansion of these 
modes in order to analyze the RC in basis of of these basis modes. This helps analytical understanding of the spiral 
galaxies more transparently. The properties of each modes can be easily understood because they are integrable. 
The corresponding coefficients C n i and a-j will determine the contributions of each modes to any galaxies. One will 
then be able to construct tables for spiral galaxies with the corresponding coefficients of each modes. Hopefully, this 
expansion method will provide us a new way to look at the major dynamics of the spiral galaxies. Analytic approach 
to the dynamics of highly flattened galaxies 

In summary, the compact expressions (trust and drum have been shown in this paper to be useful in the estimate 
of the mass density g(r < R) and g(r < R) from the exterior data at r' > R. Explicit models are presented in this 
paper. 

One also presents a more detailed derivation involving the definition of the elliptic functions E and K . Various 
useful formulae are also presented for heuristic purpose. 

One also focuses on the application in the case of modified Newtonian dynamics for two different successful models: 
the Milgrom model and the FB model. The theory of MOND appears to be a very successful model representing 
possible alternative to the dark matter approach. Nonetheless, MOND could also be the collective effect of some 
quantum fields under active investigations |lCl I22H . The method and examples shown in this paper should be of help 
in resolving the quested puzzle. 
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